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Abstract 

Cosmology is presently facing the deep mystery of the origin of the observed 
accelerated expansion of the Universe. Be it a cosmological constant, a homo- 
geneous scalar field, or a more complex inhomogeneous field possibly inducing 
effective modifications of the laws of gravity, such elusive physical entity is 
indicated with the general term of "Dark Energy" . The growing role played 
by numerical N-body simulations in cosmological studies as a fundamental 
connection between theoretical modeling and direct observations has led to 
impressive advancements also in the development and application of specific 
algorithms designed to probe a wide range of Dark Energy scenarios. Over 
the last decade, a large number of independent and complementary investiga- 
tions have been carried out in the field of Dark Energy N-body simulations, 
starting from the simplest case of homogeneous Dark Energy models up to 
the recent development of highly sophisticated iterative solvers for a variety 
of Modified Gravity theories. In this Review - which is meant to be comple- 



mentary to the general Review by Kuhlen et al. published in this Volume - I 



will discuss the range of scenarios for the cosmic acceleration that have been 
successfully investigated by means of dedicated N-body simulations, and I 
will provide a broad summary of the main results that have been obtained 
in this rather new research field. I will focus the discussion on a few selected 
studies that have led to particularly significant advancements in the field, 
and I will provide a comprehensive list of references for a larger number of 
related works. Due to the vastness of the topic, the discussion will not enter 
into the finest details of the different implementations and will mainly focus 
on the outcomes of the various simulations studies. Although quite recent, 
the field of Dark Energy simulations has witnessed huge developments in the 
last few years, and presently stands as a reliable approach to the investigation 
of the fundamental nature of Dark Energy. 
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1. Introduction 



After centuries of philosophical speculation about the origin and the phys- 
ical properties of the Universe, at the beginning of the last century cosmology 
was finally allowed to become a proper scientific discipline with the devel- 



opment of Einstein's theory of General Relativity in 1915 (Einstein, 1915) 



and with the subsequent derivation of cosmological solutions to Einstein's 



field equations by Friedmann in 1922 (Friedman, 1922). Less than a hundred 



years later, we are now provided with a well-established framework to study 
the properties of the Universe as a whole and to interpret an ever increasing 
amount of high-quality observational data that allow to continuously improve 
the constraints on a few basic parameters that fully characterize our present 
standard cosmological model. 

The latter is based on the assumption of homogeneity and isotropy of 
space encoded by the Copernican principle, and on the observation of the 
cosmic expansion that was first detected by Slipher and Hubble in the end 

This fundamental observation, 



Slipher 


1927 


Hubble 


1929 



which clearly indicated a time evolution of the Universe and posed the basis 
for the development of the Hot Big Bang cosmological scenario, removed any 
motivation for the quest of static solutions to the field equations of General 



Relativity, and led Einstein to reject his own hypothesis (Einstein, 1917) 
of a cosmological constant term that could prevent a static Universe from 
collapsing under its own self-gravity. 

The idea of a cosmological constant A acting as a sort of "repulsive force" 
and capable to counteract the attractive pull of gravity was then disregarded 
for most of the century, until new observations of galaxy correlations at large 



scales (Maddox et al. , 1990) started to indicate a tension with the predic- 
tions of a fiat matter-dominated Universe. Finally, at the very end of the 20th 
century, the extraordinary discovery that the cosmic expansion is presently 



accelerating (Riess et al. 1998 Perlmutter et al. 1999 Schmidt et al. 1998) 



suddenly revived the interest in the cosmological constant as the simplest 
possible explanation for such new observational evidence. Together with the 
wide range of astrophysical data supporting the existence of Cold Dark Mat- 



ter (CDM) as the main fraction of the total cosmic mass (see e.g. Bertone 



Hooper, and Silk, 2005 Bergstrom 2012), the discovery of the accelerated 
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expansion represents one of the observational pillars on which the presently 
accepted standard cosmological model is founded. 



Despite the remarkable success of the simple original idea of a cosmologi- 
cal constant in describing the observed properties of the accelerating Universe 
- as a consequence of which the standard model takes the name of "ACDM" 
cosmology - the theoretical roots of such idea are yet poorly defined and 
difficult to accommodate in the context of General Relativity and Quantum 



Field Theory (see e.g. Weinberg, 1989). As a matter of fact, the cosmo- 
logical constant A has to be highly fine-tuned with respect to the natural 
energy scales of the early Universe in order to provide the excellent fit to 
cosmological observations that presently still supports its success. For this 
reason, alternative explanations for the observed cosmic acceleration have 
been proposed, and are generically indicated with the term "Dark Energy" . 

Dark Energy (DE) is then simply a label with which cosmologists indi- 
cate any physical mechanism capable to provide an acceleration of the cosmic 
expansion compatible with our present observational constraints. Such pos- 
sible mechanisms - which include the cosmological constant as the simplest 
option - encompass a wide range of other alternative and more sophisti- 
cated possibilities. These include, among others, new fields and interactions 
in the Universe, cosmological models with extra dimensions, modifications of 
General Relativity, local deviations from the Copernican principle, and back- 
reaction effects of the formation of cosmic structures on the overall cosmic 



expansion (for a general and recent review, see e.g. Amendola et al. , 2012) 



Most of the present efforts of theoretical and observational cosmologists 
are devoted to the investigation of the DE phenomenon, with the aim to 
restrict the range of potentially viable scenarios for the cosmic acceleration 
and to constrain their specific parameters. In such context, several ambitious 
observational initiatives have been put in place worldwide to probe the na- 
ture of DE, and will provide complementary data of unprecedented quality 
over the next decade. These include e.g. the Dark Energy Survey (DES, 



(HETDEX, 



Abbott et al. , 2005), the Hobby- Eberly Telescope Dark Energy Experiment 



Hill et al. , 2008p , the Large Synoptic Survey Telescope (LSST, 



Ivezic et al. , 2008 ) and the recently selected European Space Agency sateUite 
mission Euclid (Laureijs et al. , 2011) that will be launched in 2020. Such 



large amount of data will have to be confronted with a wide variety of the- 
oretical proposals of ever increasing complexity and sophistication (see e.g. 
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the recent and comprehensive review of the Euchd collaboration, Amendola 



et al. , 2012) with the aim to detect possible specific observational footprints 
identifying a particular DE candidate. As a matter of fact, the detection of 
any deviation from the expected behavior of a cosmological constant would 
represent a breakthrough in our understanding of the Universe and would 
open the way for the discovery of new physics. 

The comparison between observational data and theoretical models of 
the Universe is however not a straightforward process. Besides the ever more 
complex procedures required to reduce raw data, quantify systematic errors, 
and extract meaningful cosmological information from direct observations, 
one also needs to take into account the corresponding difficulty of provid- 
ing reliable theoretical predictions for the same observable quantities. In 
fact, these often require to model highly nonlinear processes and involve the 
superposition of different physical mechanisms with potentially degenerate 
effects. 

In this respect, the use of numerical simulations to investigate the evo- 
lution of the Universe and the formation of cosmic structures beyond the 
linear regime that is readily accessible to analytical computations has proven 
to be an extremely valuable tool for the development of our understanding 
of the Cosmos. This is already true for the simplest standard ACDM model, 
but it becomes even more relevant for more complex DE scenarios for which 
one aims at identifying small deviations from the standard predictions and 
looking for such small deviations in the data. Significant progress has been 
made in the field of cosmological numerical simulations over the last decades, 
both due to the increase of the available computational power and to the de- 
velopment of efficient and sophisticated algorithms. These have allowed to 
study in detail the nature of Dark Matter and its role in driving the growth 
of cosmic structures starting from the tiny density fiuctuations generated in 
the early Universe by the infiationary accelerated expansion, and to establish 
the CDM paradigm as the main framework for the formation of galaxies and 



galaxy clusters (see e.g. the general Review by Kuhlen et al. , 2012, included 
in the present Volume). More recently, and in particular after the discov- 
ery of the cosmic acceleration, numerical simulations have also been used to 
test the nature of DE, by employing ever more sophisticated implementa- 
tions capable of capturing the characteristic features of several different and 
competing DE candidate models. Although this is a quite new and rapidly 
developing field, numerical simulations of DE scenarios beyond the cosmolog- 
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ical constant have now made sufficient progress to deserve full consideration 
as a robust and reliable approach to the investigation of the DE phenomenon. 
Therefore, cosmological N-body simulations now stand as an essential link 
between theoretical modeling and direct observations for any present and fu- 
ture collaborative initiative aimed at the study of the accelerated expansion 
of the Universe. 



The present Review is meant to provide a broad overview on the de- 
velopments and the results achieved in the field of numerical simulations for 
different DE models. The focus will be more concentrated on the conclusions 
reached by different simulation codes rather than on their numerical imple- 
mentation details. Also, due to the vastness of the topic, it will be clearly 
impossible to discuss most of the results presented in this work in full detail, 
and consequently this Review should be mainly taken as a general reference 
to address potentially interested readers to the relevant literature. A more 
general Review on cosmological N-body simulations mostly focused on the 



study of Dark Matter properties has been recently compiled by Kuhlen, Vo- 



gelsberger, and Angulo| and can be found as a separate contribution to this 



Special Issue (Kuhlen et al. 2012). 



The paper is organized as follows. In Section |2] I will present some his- 
torical outline on the role played by cosmological N-body simulations in the 
investigation of the DE phenomenon; in Section[3]l will briefly summarize the 
main classes of DE models alternative to the standard cosmological constant; 
in Section |4] I will review recent results of N-body simulations for DE mod- 
els that only modify the background expansion history of the Universe with 
respect to ACDM; in Section |5] I will then review the results of simulations 
for models where the DE also directly alters the growth of cosmic structures 
due to its density perturbations or interactions. Finally, in Section [7] I will 
provide a summary and drive my conclusions. 



2. Dark Energy and numerical simulations: some historical re- 
marks 

Numerical N-body simulations have been very successfully employed over 
the last fifty years to study the properties and the formation processes of col- 
lapsed systems in the Universe, and significantly contributed to establish the 
Cold Dark Matter (CDM) paradigm as the standard scenario for structure 
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formation (see e.g 


Aarseth 


1963; 


Peebles 


1970 


White, 


1976 


Frenk et al. 


1983 


Davis et al. 


1985; White et al 


1987; 


Navarro et al. , 


1996, 


1997 Klypin 


et al. 


19991 Moore et al. 1 


999 ISpringel et al. 


2005; Springel et al. |2008l 


Angulo et al. 201 


2). However, cosmological simulations have also played a 



major role in the discovery and in the subsequent investigation of the DE 
phenomenon. In fact, despite the undoubtable importance of the direct de- 
tection of the cosmic acceleration by Perlmutter, Riess and Schmidt (recently 
recognized also by the award of the Physics Nobel Prize 2011) it is worth to 
remind that the first observational claim of a DE-dominated universe came 
about ten years before from the comparison of the large-scale correlation of 
galaxies in the APM galaxy survey with the predictions of N-body simula- 



tions (Maddox et al. 



1990 



Efstathiou et al. 1990). 



In particular, Maddox et al. compared the correlation function extracted 



from the simulations of a CDM dominated Universe performed by [White 
et al. (1987) with the APM observational correlation function, and found 
a stark discrepancy between the two for large correlation angles, with the 
latter showing a higher level of clustering at large scales as compared to 
the numerical predictions. Shortly after, Efstathiou et al. (1990) showed 
that such large discrepancy was removed when comparing the data with 
simulations of a fiat low-density Universe with Qm ~ 0.2, where the missing 
energy for closure was given by a Cosmological Constant A. Therefore, it 
seems not inappropriate to state that the first observational evidence of a DE- 
dominated Universe was actually derived from the outcomes of cosmological 
N-body simulations. 

The connection between N-body simulations and DE investigations is 
then definitely not a new research field, although as a matter of fact it was 
only relatively recently that simulation codes suitable to explore a signifi- 
cant range of DE scenarios beyond the standard ACDM cosmological model 
started to be developed and applied. For long time, in fact, most of the 
efforts in numerical cosmology have been devoted to improve the efficiency 
and the scalability of standard N-body algorithms for the ACDM scenario. 
Such efforts have been mainly driven by the aim to reach higher and higher 
levels of detail in the description of the properties of nonlinear structure for- 
mation, as well as to include in the integration scheme the effects of baryonic 
physics (see ( 



Teyssier 


2002 


Springel and Hernquist 


2002 


Duffy et al. 



2010) and a wide range of astrophysical processes such as gas cooling, star 



formation, and feedback mechanisms from supernovae explosions and AGN 



Springel and Hernquist 


2003b|a 


Kay et al. , 


2002 


Schaye 



20041 ISijacki et all [20071 IDalla Vecchia and Schayel 120081). Alternatively, 



large N-body simulations of the standard ACDM scenario have also been 
used to develop and calibrate semi-analytic methods to populate simulated 



CDM halo catalogs with realistic galaxy samples (White and Frenk 1991| 


Lacey and Cole 1993 Kauffmann et al.[ 1993 Cole et al. 


1994| I 


CaufFmann 


et al. 


1999[ 


Somerville and Primack 1999| Springel et al. 


2001a 


De Lucia 


et al. 


2006) 





Both these approaches have driven spectacular progress in the under- 
standing of galaxy formation and evolution as well as in the capability of 
directly relating the outcomes of large numerical simulations to real obser- 
vations of galaxy and cluster populations. We are then now provided with 
a sophisticated and robust numerical machinery for simulating the evolution 
of primordial density perturbations into a wide variety of possible observable 
quantities. Nonetheless, certainly due to the excellent fit that a simple cos- 
mological constant provides to most presently available data, all such devel- 
opments have been pursued assuming a ACDM cosmology as the framework 
within which complex astrophysical processes should take place. However, 
from a theoretical perspective the cosmological constant does not appear as 
a satisfactory explanation of the DE phenomenon, and a wide range of al- 
ternative scenarios have been proposed insofar, as already briefly mentioned 
above. The attempt to include such alternative scenarios into the capabilities 
of N-body algorithms - with the aim to investigate their effects on structure 
formation processes - comes then as a natural further step in the connection 
between theoretical and observational cosmology. 



3. Dark Energy models 

It is not surprising that the discovery of the accelerated expansion of the 
Universe triggered a great deal of theoretical attempts to provide a sensible 
explanation (possibly with a lower degree of fine-tuning than the cosmologi- 
cal constant) to this mysterious phenomenon. New DE models are proposed 
almost on a daily basis (since June 1998, the number of papers containing the 
term "dark energy" in the title is about 3 thousand, corresponding to more 
than a paper every second dajQ, and often do not differ sufficiently from 
each other in their observational predictions to be possibly distinguished by 



data from www.arXiv.org 
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presently available data. A complementary approach to the development of 
specific DE scenarios based on different assumptions or on additional physical 
degrees of freedom with respect to the standard model, is that of parameteriz- 
ing our ignorance about the fundamental nature of DE with a few parameters 
quantifying possible deviations from the ACDM behavior. In both cases, in 
order to obtain realistic predictions for observables that involve - directly or 
indirectly - the nonlinear evolution of cosmic structures, it is necessary to 
include the characteristic features of each specific model or parameterization 
into the algorithms of cosmological N-body solvers. 



The range of available models and parameterizations is indeed quite large, 
including violations of large-scale homogeneity and isotropy, new dynamical 
fields, effective or fundamental modifications of the laws of gravity, and ex- 
tra dimensions. It clearly goes beyond the scope of the present Review to 
present and discuss in detail the main features of all these different exten- 
sions of the standard model, for which I refer to some specific recent reviews 



(see e.g. 


Copeland et al. 


Amendola et al. 




2012 


). ] 



2006; 


Tsujikawa 


2010 


Sapone 


2010; 


Kunz 


2012 



a sensible classification of DE scenarios should be based on the way in which 
different models can possibly affect the processes of structure formation, and 
in particular on how they are expected to modify the nonlinear collapse of 
gravitationally bound systems. Following this general principle, we can define 
three main categories of DE models: Homogeneous DE fields, Inhomo- 
geneous DE fields, and Large- Void inhomogeneous cosmologies. Far 
from trying to be complete, I will briefiy summarize the main features and 
the most common examples of these three different classes in the remainder 
of this section. 



3.1. Dark Energy as a homogeneous field 

For a homogeneous and isotropic Universe described by a Friedmann- 
Lemaitre-Robertson- Walker (FLRW) metric: 

ds"^ = -c^dt^ + a{t) [5ijdx'dx^] (1) 

where the time dependence of the line element is confined in the scale factor 
a(t), the background evolution of the Universe is encoded by the Hubble 
function H{a) = a/a which describes how the expansion rate changes as a 
function of time. Here an overdot represents a derivative with respect to 
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the cosmic time t and I assume the scale factor a to be normahzed at unity 
today. The Hubble function is then related to the relative abundance of the 
different constituents of the Universe through the Friedmann equation: 

= Qua-^ + Qra-' + QKa-' + l^pEexp {-3 ^ ^ J^"'^ cia^| , (2) 

where Qi is the energy density of the i-th component of the Universe at the 
present time in units of the critical density pcrit = SHq/ShG, and the differ- 
ent components considered are matter (M), radiation (r), curvature (K) and 
Dark Energy (DE). The equation of state parameter w{a) quantifies the ratio 
between pressure and energy densities of the DE component, and is allowed 
to be time-dependent. As one can see from Eq. [2| a cosmological constant 
corresponds to a constant value of tf = — 1, which implies a constant energy 
density of DE throughout the whole expansion history of the Universe. On 
the other hand, different constant or time-dependent values of the equation 
of state parameter would imply some evolution of the DE density and would 
consequently affect the expansion rate H{a). 

In the late Universe (i.e. sufficiently after matter-radiation equivalence at 
Zeq ~ 3x10^) the growth of linear density perturbations at sub-horizon scales 
is described, in the Newtonian gauge and in Fourier space, by the following 
evolution equation: 

6m + 2H6m = 4:T!-G (pm^m + Pde^de) , (3) 

where 5m,de = 5pm,de/pm,de is the density contrast of the matter and DE 
components. If one assumes that the DE field does not appreciably cluster at 
sub-horizon scales, i.e. if the DE component is homogeneous over the whole 
causally connected Universe, the second term on the right-hand side of Eq. [3] 
vanishes at all times since 5de = 0, and the only impact that DE can have on 
structure formation processes comes through the Hubble friction term 2H6u 
appearing on the right-hand side of Eq. [3j Therefore, a non-standard yet ho- 
mogeneous DE component characterized by an equation of state parameter 
w —1 will affect the evolution of density perturbations only in an indirect 
way through a different expansion history. Nevertheless, the impact of this 
class of scenarios on the linear and nonlinear evolution of structures can still 
be substantial as the gravitational collapse of density perturbations will oc- 
cur at different epochs depending on the evolution of the linear growth factor. 



9 



The homogeneity of the DE field at sub-horizon scales can either be taken 
as an assumption for a wide range of phenomenological parameterizations of 
the DE background evolution, or can arise as an intrinsic feature of DE sce- 
narios based on the dynamical evolution of a light scalar field (f) as in the 
case of Quintessence (Wetterich, 1988; Ratra and Peebles, 1988), k-essence 



(Armendariz-Picon et al. , 2001), Phantom (Caldwell, 2002) and Quintom 



(Feng et al. , 2005) DE models. The latter are generally characterized by a 
scalar field sound speed equal or comparable to the speed of light c, thereby 
suppressing perturbations of the DE density within the cosmic horizon, while 
DE perturbations remain frozen to a constant amplitude at super-horizon 
scales. This implies that density fiuctuations in the DE field are in any case 
present at scales comparable to the cosmological horizon even for scalar field 
models with a high sound speed ~ c (see e.g. |Ma et al.l |l999l |Bean and 



~ c (see e.g. 


Ma et al. 


Bartolo et al. 


2004 


). 1 



Dorej |2004[ [Weller and Lewisj |2003[ [Bartolo et al.| [2004j). In particular, DE 
perturbations might still change the large-scale shape of the matter power 
spectrum, thereby affecting the initial conditions for structure formation (see 
e.g. 



Ma et al. , 


1999 


Alimi et al. 


2010 



scale DE perturbations on the nonlinear evolution of structures for this class 
of models is rather small and can be expected to play a significant role only 
for very large cosmological simulations with a comoving size comparable to 
the cosmic horizon. For this reason, assuming the homogeneity of the DE 
field for this class of scenarios represents a valid approximation for a wide 
range of numerical setups, while only the recent development of extremely 
large N-body simulations for DE cosmologies (see e.g. Alimi et al. , 2010 



Rasera et al. 2010; Alimi et al. , 2012) has required to carefully take into 



account the presence of DE perturbations in the initial conditions. 



As a general example of scalar field DE models, the dynamic equation 
of a Quintessence scalar field is described by a homogeneous Klein-Gordon 
equation 

where V{(j)) is a self-interaction potential and where the DE density is given 
by Pde = 0^/2 + V{(j)). Different choices of the function V{(f)) will then 
determine different evolutions of the DE density and will affect in specific 
ways the expansion history and consequently the growth of cosmic structures. 
Some of the most widely used forms of the function V{(j)) include runaway 
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potentials as e.g. an inverse power-law (Ratra and Peebles 



1988) 



(5) 



or an exponential (Lucchin and Matarrese, 1985 Wetterich, 1988): 



1^(0) = Ae-°^ , 



(6) 



as well as con/inmg' functions as e.g. the SUGRA potential (Brax and Martin 



1999) arising naturally within supersymmetric theories of gravity: 

Vi<l>) 



(7) 



In Eqs. |5]|7] the scalar field 
Planck mass Mpi 



has been expressed in units of the reduced 
l/\/87rG and is therefore dimensionless. These three 



potentials represent the most widely used choices for Quintessence models as 
they provide viable expansion histories and scaling solutions that make the 
cosmological evolution largely independent from the scalar field initial condi- 
tions (see e.g. Ferreira and Joyce, 1998), and have been widely investigated 



through N-body simulations. 

Cosmological models characterized by a vector field, rather than a scalar. 



playing the role of DE have also been recently proposed (Beltran Jimenez and 



Maroto, 2008). In such scenarios, cosmic acceleration is driven by the kinetic 



energy of the vector field, without resorting on any arbitrary choice of a po- 
tential function. Despite the vector nature of the DE field, the energy density 
of its spatial components dilutes faster than matter with the cosmic expan- 
sion, and is therefore negligible for the evolution of the late Universe. These 
models therefore behave similarly to scalar field DE cosmologies, inducing a 
modified background expansion history without significant sub-horizon per- 
turbations of the DE density, although the fundamental mechanism behind 
the accelerated expansion is different from standard Quintessence scenarios. 

A different possibility, already mentioned above, is to assume a priori the 
homogeneity of the DE field and describe its time evolution by phenomenolog- 
ical parameterizations of the DE equation of state parameter w{a). Several 
different options have ben proposed in the last years, either based on the 
behavior of wia) at low redshifts, as for the case of the so-called Chevallier- 



Polarski-Linder (CPL, Chevallier and Polarski 2001 Linder 2003) parame- 
terization 

w{a) = wq + Wa{l - a) (8) 
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where Wq and Wa are constants, or assuming as main parameters the relative 
abundance of DE at the present epoch (1 — f^n) and at early times (I^ede), as 
for the case of the Early Dark Energy parameterization of Wetterich] ( 2004[ ): 



w(a) = , , , where b = — 7—^ — . ^ • (9) 

l + 61n(l/a) In + In i-f^M 



All these different scenarios and parameterizations significantly affect the 
growth of density perturbations both in the linear and nonlinear regimes, 
and have been extensively investigated with N-body simulations over the 
last decade, as will be discussed in Section |4j 

3.2. Dark Energy as an inhomogeneous field 

If DE is associated to some new physical degree of freedom rather than 
to a cosmological constant, it is natural to consider also its spatial fluctu- 
ations and its possible interactions with other components of the Universe. 
The assumption of homogeneity discussed above might therefore be a rea- 
sonable approximation at sufficiently small scales for a wide range of DE 
scenarios characterized by a large sound speed of the DE fluid and by the 
absence of substantial direct interactions of DE besides gravity, but certainly 
fails in describing the most general possible case of a DE field beyond A. A 
large variety of DE models featuring significant perturbations at sub-horizon 
scales and/or substantial interactions with matter or gravity have been pro- 
posed in the last years, and generically form the class of inhomogeneous DE 
cosmologies. 

A first example of such models is given by the Clustering DE scenario 
(e.g. Creminelli et al. 2009 2010 Sefusatti and Vernizzi, 2011), where a 



general k-essence scalar field is simultaneously characterized by an equation 
of state parameter generally different from —1 and by a "cold" sound speed 
Cg ~ 0. As a consequence, DE can cluster also below the horizon and source 
gravitational potentials at scales relevant for the formation of cosmic struc- 
tures. This corresponds to the case 5de 7^ that was discarded above under 
the assumption of homogeneity, which implies that the net potential for the 
growth of CDM density perturbations will include also the contribution of 
DE perturbations, according to the full form of Eq. [3} for which one can 
write: 

5m + 2HSm = 47rGpM ( Sm + ^^^^de) • (10) 
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From Eq. 10 one can clearly see that DE perturbations will substantially 
affect the evolution of CDM structures only at late times, when the back- 
ground DE density becomes important as compared to the CDM one. Also, 
it is interesting to notice how an observer ignoring the clustering properties 



of DE could interpret the evolution of perturbations determined by Eq. 10 
as a modification of gravity emerging in the late Universe. The example of 
Clustering DE models then already clearly shows how a fundamental dis- 
tinction between a DE degree of freedom and a modification of the laws of 
gravity at astrophysical scales results impossible whenever one allows for spa- 
tial perturbations in the DE field: the specific clustering properties of a DE 
field can in general mimic deviations from the expected behavior of standard 
gravitational instability processes induced by a modified force law. 

The fundamental degeneracy between these two different perspectives be- 
comes even more evident for the case of DE fields featuring direct interac- 
tions with matter, for which a formal correspondence to modified theories 
of gravity through a conformal transformation of the metric can be explic- 



itly demonstrated (see e.g. Pettorino and Baccigalupi, 2008). Interacting 



DE models and Modified Gravity theories therefore represent a unique class 
of cosmological scenarios beyond ACDM for which structure formation pro- 
cesses are in principle modified both by a non-standard evolution of the 
background expansion history and by the specific clustering and interaction 
properties of the new degrees of freedom associated to the DE sector of the 
Universe. Such cosmologies are in fact generically characterized by the exis- 
tence of fifth-forces mediated by these new degrees of freedom, whose spatial 
range and universality depend on the specific model under consideration. A 
detailed overview and classification of Interacting DE and Modified Grav- 
ity models goes beyond the scope of the present Review, and I refer the 
interested reader to some excellent recent publications which provide a self- 



consistent and comprehensive overview on these scenarios (see e.g. Tsujikawa 



2010 De Felice and Tsujikawa 2010 Amendola et al. , 2012 and references 



therein). For what concerns the aims of this work, it is sufficient to identify 
the few main features that determine how different specific models belonging 
to this class of cosmologies directly affect the growth of density perturbations 
in the linear and nonlinear regimes. 

As already mentioned above, a general feature of Interacting DE and 
Modified Gravity models is the existence of a fifth-force of nature, mediated 
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by the scalar degree of freedom associated to DE. In the most general case, 
such fifth-force can be described as an additional term in the acceleration 
equation of a massive test particle representing a fluid element of a given 
cosmic component i: 

= -V<f-A(0)W(/. (11) 

where the standard gravitational potential $ is determined by the Poisson 
equation: 

V^«l> = 47rGj]pA (12) 

j 

with j ranging over all the different clustering components of the Universe. 
The additional scalar potential S(j) obeys a modified non-linear Poisson equa- 
tion of the form: 



^^Scj) = F{5(t)) + 87iG(3j 



(13) 



with F a generic function of the scalar field perturbation S(j). As one can 



see from Eqs. 11 13 the choice of the coupling functions A(0) and the form 
of the function F{6(j)) will determine the configuration of the scalar pertur- 
bations 5(j) and the related fifth-force experienced by massive particles. The 



formulation presented above and described by Eqs. 11 13 is rather general, 
and covers a wide range of different models of Interacting DE and Modified 
Gravity. 



As a first main classification of such scenarios, one can then start distin- 
guishing between models featuring a universal coupling (i.e. /3i((/>) = /3(0) Wi) 
and models with species- dependent couplings. The former case, correspond- 



ing to modified gravity theories as e.g. f{R) gravity (see e.g. Hu and Saw- 



icki, 2007 De Felice and Tsujikawa 2010 and references therein). Extended 



Quintessence models (Baccigalupi et al. 2000 Perrotta et al. 2000 Pettorino 



et al. , 2005 Pettorino and Baccigalupi 2008), higher-dimensional theories of 



gravity as e.g. DGP (Dvali et al. 



2000), or the recently proposed Galilean 



(Nicolis et al. , 2009), Symmetron ( Hinterbichler and Khoury, 2010 Hinter 



bichler et al. , 2011) and Dilaton (Gasperini et al. , 2002) cosmologies, requires 



that the fifth-force be suppressed in high-density environments in order to 
evade solar system constraints on possible deviations from General Relativity 
(see e.g. Bertotti et al. 2003; Will 2005). This suppression can be realized 



with a variety of screening mechanisms, as e.g. the Chameleon (Khoury and 
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Weltmanl 120041), the Vainshtein (IVainshteinl fT972j IDeffayet et all 120021) or 



the Symmetron (Hinterbichler and Khoury, 2010), which all rely on relatively 



large fluctuations {\5(j)\ ~ 0{1) or larger) of the scalar field (or of its deriva- 
tives) between high-density regions and the cosmic low-density environment. 
Such large perturbations can arise e.g. when nonlinearities are present in the 



function F[5(j)) appearing in Eq. 13, which then requires quite sophisticated 



algorithms to be properly solved for an arbitrary matter distribution (5m (^ , x) 
within newtonian N-body codes, as will be discussed in Section |5| This is for 
instance the case of f{R) theories of gravity in the Hu and Sawicki parame- 
terization, for which 6(j) = fji^ df{R)/dR and F^Scf)) = F{fji) oc -R(/r) — R 
with the relation between and R given by: 



„ ci 



n+1 



(14) 



with n, Ci and C2 constants. 

On the other hand, if one allows for non-universal couplings (as first 
proposed by Damour et al. , 1990), solar system constraints can be easily 



evaded without resorting on any screening mechanism by simply assuming 
the coupling to baryons (3b{4>) to be highly suppressed. This second option 
corresponds to the general class of Coupled DE models where a non- vanishing 



coupling to CDM particles (Wetterich, 1995 Amendola 2000, 2004) or to 



massive neutrinos (as for the Growing Neutrino scenario, Amendola et al 



2008 ) provides viable cosmological expansion histories and a possible solution 



to the fine-tuning problems of the cosmological constant. For this class of 
models, the function -F(50) in Eq. [Tsjis related to the derivative of the scalar 
self-interaction potential dV/ d(j), and for sufficiently flat potentials (which are 
anyway required in order to provide an accelerated expansion of the Universe) 
can be safely discarded compared to the matter density perturbations, such 



that Eq. 13 reduces to: 



V^(50fii^87rG/3j((/))5j 



(15) 



and for the case where only one species is coupled to DE one gets: 
V^6<j) ^ 87rGA(0)(5, = 2A(0)V'<I', 



(16) 



where <I>j is the standard gravitational potential generated by the coupled 
matter component i. From the previous equation, one immediately gets that 
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(50 ~ 2/3j(0)$j and therefore from Eq. 11 the fifth-force acting on a coupled 
particle will be simply proportional to standard gravity by a factor 2l3f. More 



general scenarios featuring a coupling with two (Baldi, 2012b) or multiple 



(Brookfield et al. , 2008) CDM fluids have also been recently proposed, for 
which the previous arguments apply separately to the fifth-force generated 
by each individual coupled component. 

For the case of a non-universal interaction between DE and other fluids 



as a conse- 



in the Universe, an additional acceleration term appears in Eq. VI 
quence of momentum conservation in the coordinate frame of the minimally 
coupled species (i.e. those species for which the couphng to DE vanishes, see 
e.g. 



Amendola 2000 


Maccio et al. 


2004 


Baldi et al. 


2010 



term is in general proportional to the velocity vector of a test particle and has 
been therefore termed "friction" or "drag" term in the literature. The full 
acceleration equation of a coupled particle in the Einstein frame for Coupled 
DE models with a non-universal coupling then reads: 



Vi 



- V$- 2/3,(0) 5^/3,(0) V% 



(17) 



which for a self-consistent N-body implementation requires to separately 
solve for the gravitational potential of each differently-coupled matter com- 
ponent of the Universe. It is also interesting to notice here how the sign 
of the friction term depends on the relative signs of the scalar fleld back- 
ground velocity and of the coupling function /3i(0). This peculiar form of 
the friction term can determine a quite broad phenomenology of interacting 
DE models at the level of linear and nonlinear structure formation, as will 
be discussed in Section [5] below. 

3.3. Large- Void models 

A further possibility to account for the accelerated expansion of the Uni- 
verse without invoking a Cosmological Constant (and in this speciflc case 
even without resorting on any other DE fleld) is to drop the assumption of 
large-scale homogeneity encoded in the Copernican Principle and consider 
the possibility that the observed accelerated expansion be just an apparent 
effect due to a strong local deviation from homogeneity (see e.g. Mustapha 



eFaLj ITomitaj [200T| |Wiltshirel [20071 ICarcia-BeUido and Haugboelle 



2008). In particular, an observer sitting near the center of a very large un- 



derdensity would observe an apparent acceleration of the Universe due to the 
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different expansion rate of the void at different distances from its geomet- 
rical center. This class of scenarios goes under the name of Large- Void or 
Lemaitre-Tolman-Bondi (LTB) cosmologies as they derive from the general 
spherically symmetric space-time metric first studied by Lemaitre, Tolman 
and Bondi ( |Lemaitre| p33| |Tolman| \l93^ |Bondi[ [l947| : 



-de + ^''^'''y''' + A\r, t) (de' + sin^ edcl)') , 
1 — k{r) ^ ' 



(18) 



where A and A! are functions of time and of the radial coordinate from the 
center of symmetry of the system. 

Although LTB models require a very large size of the density void (~ 2 
Gpc or larger) in order to possibly explain the observed accelerated expan- 
sion of the Universe without resorting to any additional DE field, they have 
attracted significant interest in the last years due to their simplicity and to 
the wide range of possible observational features that they provide and that 
could become directly observable with the next generation of surveys (see 



e.g. Quercellini et al. 2010). 



Viable large- void LTB cosmologies can be described by a four-parameters 
model of the void density profile ^y\{j) and of the radial Hubble rate H^ij) 
according to the equations 



see 



Garcia-Bellido and Haugboelle 20081 for 



more details): 



l]M(r) = 1 + (fiin - 1) 



1 — tanh \{r — rg) /2r] 
1 + tanh [ro/2Ar] 



^M(r) /^^K(r) 



(19) 
(20) 



where the four free parameters are the overall expansion rate i^o, the un- 
derdensity at the center of the void f2in, the radius of the void ro, and the 
transition width of the void profile Ar which defines how the profile matches 
from the inner value f2in and the density parameter at infinity which is as- 
sumed to be f2M(cxD) = 1. Such class of models affects the growth of density 
perturbations due to the space-dependence of the cosmic density '^y\{r^ which 
for very large voids will still be approximately constant over the scales of den- 
sity perturbations collapsing into bound structures before the present epoch, 
but will significantly vary over different regions of the presently observable 
Universe. Despite large-void LTB models have recently started to show some 
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tension with geometric probes of the expansion history of the Universe (Zu- 



malacarregui et al. , 2012), the study of their effects on structure formation 



processes with the aim to identify possible observational footprints of a large 
void in the statistical properties of large-scale structures has recently at- 
tracted significant interest, and will be briefly discussed in Section |6l 



4. Simulating a Dark Energy background expansion 

I now move to discuss how the different cosmological scenarios beyond 
ACDM that were introduced in Section [3] have been investigated by means of 
dedicated N-body simulations for what concerns their effects on the forma- 
tion of nonlinear cosmic structures. I start such review from the simplest case 
of homogeneous DE models for which, as I explained above, the only effect 
on the growth of density perturbations comes from a modified background 
expansion that changes the linear growth factor through the Hubble friction 
term of Eq. [3| unless the simulated volume is so large to require a proper 
sampling of DE perturbations at scales comparable to the cosmic horizon. 
Consequently, cosmological simulations aiming at studying the evolution of 
structures in the context of these scenarios need to implement in their nu- 
merical algorithms only a proper modification of the expansion history H{z). 



The first simulations of homogeneous DE models with a constant equa- 
tion of state parameter w ^ —1 have been performed by Ma et al. (1999) 
using a Particle-Particle/Particle- Mesh code to evolve 128"* particles within 



a periodic cosmological box of 100 Mpc aside. The work of |Ma et al. fo- 



cuses mainly on the detailed shape of the nonlinear matter power spectrum 
in constant-w DE models, providing a fitting formula based on the specific 
growth factor of the different DE cosmologies. Soon after. Bode et al. (2001) 
performed a large suite of N-body simulations with 512^ particles in a 1 Gpc 
periodic box for a variety of cosmologies, including also one DE model with 
w = —2/3, and investigated the evolution of the cluster mass function in the 
different scenarios, finding that the DE model shows a slower evolution of the 
cluster abundance, thereby determining a larger number of clusters at high 
redshift when a common normalization of the linear perturbations amplitude 
at z = is assumed. 

The first simulations of homogeneous DE models with a variable equation 
of state parameter w ^ const, have been performed by |Klypin et al. (2003) 
using a modified version of the Adaptive Mesh Refinement (AMR) code ART 
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(developed by Kravtsov et al. , 1997). In their work, Klypin et al. investi- 
gated a few test models with either a constant equation of state w > —1 or 
a variable equation of state w{a) corresponding to the dynamical evolution 
of a Quintessence field with the inverse power-law and SUGRA potentials 
of Eqs. |5] and [7j For their simulations, Klypin et al. adopted a common 
normalization of the linear power spectrum of all the different cosmologies 
with the standard ACDM value of at z = 0, similarly to what previously 
done by Bode et al. As we will see later on, the choice of the linear normal- 
ization is a critical issue in the comparison of different cosmological scenarios 
with N-body simulations. The outcomes of these first runs showed that no 
significant difference was present at z = among the various models in sev- 
eral observable quantities like the nonlinear matter power spectrum P{k), 
the CDM halo mass function A^(> M), and the circular velocity function 
(i.e. the number of halos as a function of their maximum circular velocity). 
A significant scatter among the models could instead be noticed at higher 
redshifts, with the non-standard DE cosmologies systematically showing a 
higher number of halos as compared to ACDM both in the halo mass func- 
tion and in the circular velocity function, with an enhancement increasing 
with halo mass (see Fig. [T| left panel). Furthermore, the DE models did also 
show a higher amplitude of the matter power spectrum at all scales for z > 0, 
consistently with the slower growth rate induced by the background scaling 
of the DE density. 
In the same work, 



Klypin et al. also investigated the inner structure 



of CDM halos in their various DE cosmologies by means of high-resolution 
zoom re-simulations of some of the most massive halos identified in the basic 
cosmological runs. This allowed them to show that the density profiles of 
CDM halos in homogeneous DE models still follow a Navarro-Frenk- White 
(NEW, [Navarro et aL| [l997| profile 

p(r) 5* 



Peru (r/rs) ■ {l + r/rsf 



(21) 



but are systematically more concentrated than their ACDM counterparts, 
with a smaller value of the scale radius r^. According to the concusions of 



Klypin et al. this effect is likely due to the earlier formation redshift of halos 



in the DE scenarios, a picture which is again consistent with a slower growth 
rate of density perturbations for models with a common linear perturbations 
normalization at the present epoch. Such effect determines a higher normal- 
ization for the concentration-mass relation in DE cosmologies as compared 
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Figure 1: Plots from Klypin et al. (2003) - Left panel: Comparison of the HMF for a ACDM 



cosmology (red) and a Quintessence model with an inverse power-law self-interaction po- 
tential (black) as extracted from N-body simulations with a box size of 160 Mp c//t and 
a mass resolution of m = 2 x Mq/H. The other DE scenarios considered by Klypin 
lie in between these two extreme models. The mass functions are practically in- 
but at higher z the DE cosmologies show a higher number of 
Right panel: The concentration-mass relation at 



let al.l 



distinguishable at z 
massive halos as compared to ACDM 
z = computed from the same set of simulations considered in the left panel. CDM halos 
in DE cosmologies result more concentrated due to their earlier formation epoch. 
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to ACDM (see Fig. [T| right panel). Additionally, they also found that the 
DE cosmologies systematically show a higher number of CDM satellite halos 
within massive collapsed structures, and that this effect seems to correlate 
with the higher circular velocity of CDM main halos in DE models as com- 
pared to ACDM. 



A complementary approach was followed soon after by hinder and Jenkins 
( 2003[ ), who investigated DE models with a parametrized equation of state 



w{a) using the CPL parameterization of Eq. [8] with wq and Wa chosen to 
best fit Quintessence models with a SUGRA potential. With this approach, 
hinder and Jenkins] ran a series of cosmological simulations with a modi- 



fied version of the Tree-PM code GADGET ([Springel et aL| |2001b[ ) within 
a somewhat larger box size as compared to Klypin et al. (2003) in order to 
better sample the high-mass tail of the halo mass function. This study also 
adopted a common normalization of the different models to the same cxg at 
z = 0, and found consistent results with the earlier outcomes of [Klypin et al 



no significant deviations from the standard ACDM model at 2; = 0, and a 
systematically larger abundance of CDM halos - especially at large masses -~ 
for the DE cosmologies at higher redshifts as compared to ACDM. Addition- 
ally, hinder and Jenkins] showed that the the standard fitting formulae for 
the Halo Mass Function (HMF), and in particular the Jenkins et al. (2001) 
formula, still provide a good fit to the simulated HMF in DE cosmologies at 
any redshift, provided the correct growth rate of linear density perturbations 
is used in the fit. This result indicated that the universality of the HMF is 
broadly preserved in homogeneous DE models at least at the ~ 20% level, 
and that the differences in the abundance of CDM halos at high z in DE cos- 
mologies are fully captured by the different linear growth factors. A similar 
study was also performed by Lokas et al. (2004) restricting to the case of a 
constant equation of state w ^ —1, finding results in general agreement with 
these earlier claims. 

A more detailed investigation of the HMF in homogeneous DE cosmolo- 
gies has been carried out much more recently by e.g. Courtin et al. (2011) 
with higher-resolution simulations, finding evidence of deviations of the HMF 
from a universal behavior at the level of about 10% for the case of a Quintessence 
model with an inverse power-law potential. The comparison of these results 
already shows how the improvements in the simulations accuracy and dynam- 
ical range have allowed to detect progressively finer details of the imprint of 
DE on structure formation processes. 
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The effects of a homogeneous DE field on the internal properties of cluster- 



size halos was then studied in much finer details in Dolag et al. (2004) by 



running high-resolution zoom re-simulations of the 17 most massive halos 
identified in a fiducial ACDM cosmological run within a range of different ho- 
mogeneous DE cosmologies. For the DE models they considered a constant- 
w cosmology with w = —0.6 and two variable-w models corresponding to 
Quintessence scenarios with inverse power-law and SUGRA potentials, both 
with the same value of the equation of state parameter at z = 0, wq = —0.86, 
and still assuming the same cxg normalization of all the models at z = 0. With 



such setup Dolag et al. investigated the variation of the concentration param- 
eter c = Tyir/rs (where is the halo virial radius) in their high-resolution 
halo sample within the different cosmological models, finding that the over- 
concentration of halos in DE cosmologies at z = already highlighted in 



the early results of Klypin et al. (2003) can be related to the different linear 



growth factors through a simple scaling relation given by 



ACDM 




D 



DE, 



^coU) 



(22) 



where Cq is the concentration parameter at z = for a 10^^ Mq/Zi halo, 
D+{z) is the growth factor, and Zcou is the collapse redshift of the halo. The 
same study also showed that the mass dependence of the concentration-Mass 
relation is not significantly affected in homogeneous DE models, which allows 
to derive the concentration parameter at 2; = within DE cosmologies at any 
halo mass using Eq. |22] once the concentration-Mass relation is sufficiently 
tightly calibrated for the standard ACDM case. 

As a follow-up of this study, Meneghetti et al. ( 2005a|[b ) studied the strong 
lensing efficiency of the 17 clusters simulated by Dolag et al.| by means of ray- 
tracing techniques, finding that the higher concentration of clusters in the 
DE cosmologies determines a higher lensing efficiency as compared to ACDM, 
although this effect also crucially depends on the choice for the normaliza- 
tion of the linear matter power spectrum. In fact, a different normalization 
choice - assuming e.g. a common amplitude of density perturbations at the 
last-scattering surface zis ~ 1100 - would result in the opposite trend for 
all the main effects of DE cosmologies discussed so far, including a lower 
halo concentration at 2 = as compared to ACDM, and correspondingly a 
lower efficiency of clusters as strong gravitational lenses. A similar study was 
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also performed soon after by Maccio (2005) making use of the simulations 
of Klypin et al. (2003), leading to consistent results with the earlier study of 
Meneghetti et al. 



The issue of the power spectrum normalization was discussed also in 



Kuhlen et al. (2005), that investigated a series of DE models with constant 
equation of state w 7^ — 1, extending for the first time the analysis to the 
case of w < —1, generally indicated with the term "Phantom" DE. Besides 
showing that an equation of state parameter more negative than the ACDM 
value with a common normalization of the linear power spectrum to the same 
(Tg at 2; = determines the opposite trend in the resulting HMF and halo 
concentration as compared to the w > —1 case, this study also explicitly 
showed that such trends are in any case reversed if one assumes a common 
normalization of all the cosmologies to the amplitude of scalar perturbations 
at last-scattering (see Fig. |2| left). This result confirms and significantly re- 
inforces the early conclusion that the nonlinear effects of homogeneous DE 
cosmologies as compared to ACDM are mainly driven by the different evolu- 
tion with redshift of the linear perturbations amplitude in the DE cosmologies 
due to their different growth factors D+{z). 

The evolution of the baryonic component of the Universe within N-body 
simulations should be treated taking into account the coUisional nature of 
baryons as opposed to the coUisionless nature of CDM particles. There- 
fore, a variety of methods have been developed to solve the hydrodynamical 
equations of a perfect gas fluid, along with the solution of the gravitational 



interaction of masses (see e.g. Teyssier, 2002 Springel, 2010, 2011). Addi 



tionally, a number of non-adiabatic astrophysical processes can be included 
in the simulations, ranging from the radiative cooling of the gas and the 
following formation of stars, to the feedback provided by supernovae explo- 
sion and/or by the accretion of gas onto supermassive back holes, to the 
interaction between the gas and large-scale magnetic fields. The former and 
simpler approach generally goes under the name of adiabatic or non-radiative 
hydrodynamical simulations, while the latter is referred to as radiative hy- 
drodynamics. 

The first attempt to include hydrodynamical processes in cosmological 
simulations of homogeneous DE cosmologies was performed by |Maio et~al 



(2006), who studied the formation of the early gas clouds responsible for 



the reionization of the Universe in a variety of DE cosmologies, by means of 
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Figure 2: Left, plot from Kuhlen et al. (2005)- The concentration-mass relation at z = as 
extracted from a series of N-body simulations of DE cosmologies with w ^ —1, including 
also the case of "Phantom" DE, u; < — 1. The plot shows the crucial role played by the 
normalization choice of the linear perturbations amplitude: on the left, the concentration- 
mass relations for a w = —0.5 DE model (dotted) and a, w — —1.5 "Phantom" DE model 
(dashed) show that these types of DE scenarios determine respectively an increase and a 
decrease of the average halo concentration as compared to a standard ACDM cosmology 
(solid). On the right, the same plot shows the opposite trend for simulations where a 
common normalization of linear density perturbations at last scattering has been chosen. 



Right, plot from Alimi et al. 



The ratio of the nonlinear matter power spectrum 
in different DE cosmologies over the standard ACDM extracted from simulations 

with different final values of the linear perturbations amplitude. The dashed and dot- 
dashed curves represent the case of Quintessence cosmologies with an inverse power-law 
and a SUGRA potential, respectively. Different curves refer to different epochs, from top 
to bottom a = 0.3,0.5, 1. The plot shows that the maximum deviation from ACDM is 
obtained at intermediate scales k ^ 1/i/Mpc, and for low rcdshifts. 
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radiative simulations including gas cooling. In this work it was found that 
the earlier formation of structures that characterizes DE models with w > —1 
applies also to gas clouds that can then induce an earlier reionization epoch 
as compared to ACDM, a result that looked very appealing at the time due 
to the high value of the reionization redshift derived from the first-year data 



of the WMAP satellite (Spergel et al. , 2003). The same study also showed 



that a running spectral index of the primordial power spectrum ns{k) might 
significantly mitigate this effect, simply by suppressing the small-scale power 
with respect to the large-scale normalization of the linear perturbations. 

Radiative hydrodynamical cosmological simulations were also employed 
by several other authors to study the structural properties of galaxy clusters. 



Aghanim et al. (2009) performed a series of hydrodynamical runs with gas 
cooling for a range of homogeneous DE models with constant and variable w 
to study the impact of DE on the scaling relations between cluster masses and 
different mass proxies such as the cluster X-ray luminosity and the Sunyaev- 
Zeldovich (SZ) signal, finding that homogeneous DE does not significantly 
alter the standard scaling relations and concluding that the use of standard 
ACDM scaling relations also for homogeneous DE models seems generally 
appropriate. 



A similar analysis was performed by De Boni et al. (2011) 



( 


2011 


); 


De Boni 



et al. (2012), that studied the concentration-Mass relation, the luminosity- 
temperature relation, and the baryon fraction of clusters in hydrodynamical 
simulations of DE models including gas cooling and star formation, finding 
results in general agreement with previous claims. 



The impact of homogeneous DE on the nonlinear matter power spectrum 
was then investigated in detail by e.g. Ma (2007); Francis et al. (2007); 



Casarini et al. (2009). In particular both Francis et al. and Casarini et al. 
found that the nonlinear matter power spectrum of DE models with a variable 
equation of state parameter w{a) can be derived from the nonlinear power 
spectrum of constant-w models with an accuracy down to ~ 1% through a 



transformation involving only background quantities. Alimi et al. (2010) then 
also investigated the nonlinear matter power spectrum in some specific DE 
models selected to best fit background and linear perturbations observational 
data. This implies that the as normalization at z = of the different models 
be different, and generally lower, in DE cosmologies as compared to ACDM. 



With such setup [Alimi et al.| found that the maximum deviation of the DE 



power spectra with respect to ACDM occurs at intermediate scales around 
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~ 1 h/Mpc (see Fig. |2| right panel). Such behavior has been subsequently 
broadly confirmed also by Fedeli et al. (2011) for different choices of the 



homogeneous DE evolution with a similar setup, by means of cosmological 
simulations including also radiative hydrodynamical processes as gas cooling 
and star formation. [Fedeli et al. additionally showed that star formation 
efficiency is generally reduced in DE cosmologies (consistently with the earlier 



results of De Boni et al. 2011). 



This maximum deviation from the ACDM matter power spectrum at in- 
termediate scales (with an amplitude up to 40% for the most extreme model 



considered by Alimi et al. ) appears to be mostly driven by the different Ug 



normalization of the various cosmologies that was adopted both by Alimi 
et al. (2010) and Fedeli et al. (2011) (a similar feature occurs, for example, 
also for the case of Coupled DE models with high-z normalization, see e.g. 



Baldi, 2012c, and the related discussion above). In fact, although other stud- 



ies employing a common ag normalization at z = (such as e.g. Ma, 2007) 



did also find a quahtatively similar effect, its amplitude results much weaker, 
with a maximum detected deviation of the order of a few percent. However, 
such small residual deviation in the nonlinear matter power spectrum found 
for simulations with the same ag normalization represents a very important 
result as it demonstrates how the full nonlinear matter power spectrum can- 
not be uniquely determined with arbitrary precision by the amplitude and 
shape of the linear one. More specifically, this result shows that two cos- 
mological models with the same normalization of the linear matter power 
spectrum at the present epoch but with different growth histories can in 
principle be distinguished from each other through their nonlinear power 
spectra, although the deviation is expected to be small and consequently 
particularly difficult to detect. 

Another relevant effect of homogeneous DE models on observable quan- 
tities that has been investigated through N-body simulations concerns the 
Baryon Acoustic Oscillation (BAO) peak in the correlation function of col- 



lapsed halos. Jennings et al. (2009) carried out a series of large CDM-only N- 



body runs with a modified version of the TreePM code GADGET-2 (Springel 
2005[) within a box of 1500 Mpc/h aside filled with 646^ CDM particles for 



a range of homogeneous DE models with a parametrized equation of state 
w{a). For this study, Jennings et al. adopted a different parameterization 



with respect to the ones introduced in Section [3} following the conclusion 
(Bassett et al. , 2004) that the CPL parameterization does not reproduce 
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with sufficient accuracy the evolution of Quintessence models at high red- 
shifts. Instead, they employed a four-parameter parameterization proposed 



by Corasaniti and Copeland (2003) that provides a better fit to the real 



equation of state evolution for a wide range of Quintessence cosmologies. 
The outcomes of the simulations by [Jennings et al.| concerning the matter 
power spectrum and the HMF in DE models were found to be in good agree- 
ment with previous findings. Additionally, this work investigated for the first 
time the effects of homogeneous DE on the properties of BAOs, finding that 
even DE models with a significantly rapid evolution of the equation of state 
parameter at relatively low redshifts do not imprint any significant shift in 
the location of the BAO peaks as compared to ACDM, thereby making it 
difficult to detect an evolution of the DE equation of state through measure- 
ments of the BAO scale (see Fig. (sj left panel). 



The case of the Early Dark Energy (EDE) parameterization of Eq. [9 
has been treated separately from Quintessence scenarios and from the GPL 
parameterization. In two independent and almost contemporaneous works. 



Francis et al. (2008a) and Grossi and Springel (2009) investigated a range 



of EDE cosmologies by means of GDM-only simulations performed with two 
independently-developed modified versions of the N-body code GADGET, 
with a particular focus on the impact of EDE on the HMF. Both these works 
consistently found that the HMF at 2; = is only mildly affected by the 
existence of a non- vanishing fraction of DE at early times, and that the uni- 
versality of the HMF shape encoded by standard AGDM analytical formulae 
as e.g. the Jenkins and Warren ffiting functions (Jenkins et al. , 2001 War- 



ren et al. , 2006 ) is preserved in EDE cosmologies at least at the level of 



10 — 15% accuracy. These results are in contrast with previous claims by 



Bartelmann et al. (2006) based on a spherical collapse treatment of the for- 



mation of GDM halos in EDE cosmologies, which found a significant change 
in the linear overdensity at collapse 5c in the presence of an EDE compo- 



nent, and with the subsequent derivation by Fedeli and Bartelmann (2007) 



of a corresponding significant enhancement in the strong lensing efficiency 
of clusters within EDE cosmologies. Such discrepancy has been further dis- 



cussed by Francis et al. (2008b) who showed how under the assumption of 



small DE perturbations at astrophysical scales (which is the main assumption 
for homogeneous DE cosmologies and that was implicitly assumed in both 
the numerical studies mentioned above) a value of the overdensity parameter 
5r close to the standard AGDM value of 5r = 1.686 is restored. 
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Figure 3: Left, plot from Jennings et al. (2009 ): The nonlinear matter power spectrum 
divided by a smooth power spectrum (i.e. a power spectrum without baryonic wiggles) for 
the standard ACDM cosmology (black solid line) and a Quintessence model with a SUGRA 
potential (red triangles) . Although the amplitude of the baryonic acoustic oscillation signal 
is amplified in the Quintessence model, the location of the peaks coincides with the ACDM 



case (indicated by the vertical dashed lines) within a 5% accuracy. Right, plot from Grossi 
and Springel 12009 ): The velocity function A^(> a) as a function of the halo line-of-sight 



velocity dispersion a as extracted from a ACDM simulation (green line) and a series of DE 
cosmologies, including two distinct EDE models (orange and red lines). The grey-shaded 
area corresponds to the gap between the fiducial ACDM cosmology assumed in this work, 
with (Tg = 0-8, and a ACDM model with a higher cs = 0.9. 



The study of [Grossi and Springel] also investigated the concentration- mass 
relation in the context of EDE cosmologies, and the velocity function A^(> a), 
which is a conceptually similar observable to the HMF where CDM halos are 
counted based on their line-of-sight velocity dispersion a rather than by their 
total mass. Such analysis led to the interesting conclusion that the velocity 
function of EDE cosmologies at redshifts around z ~ 1.5 mimics a ACDM 
velocity function for a standard cosmology with a higher normalization of 
the linear perturbations amplitude cxg (see Fig. |3| right panel). This result 
indicates that the presence of an EDE component might be detected through 
the determination of an excessively large value of cxg from the line-of-sight 
velocity dispersion of high-z clusters. 

The case of the Vector DE models briefly mentioned in Section [3] has 



also been recently investigated with N-body simulations by Carlesi et al. 
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(2011, 2012). In their set of simulations, Carlesi et al. investigated the im- 
pact of Vector DE models on a number of observables as e.g. the cluster 
number counts as a function of redshift, the HMF, the distribution of cos- 
mic voids, and the structural properties of collapsed halos encoded by the 
concentration, spin, and shape parameters. As a main conclusion of their 



analysis, Carlesi et al.| showed that even though the large-scale properties of 



structures evolve quite differently in Vector DE cosmologies as compared to 
ACDM, such deviations are mainly driven by the different evolution of the 
background cosmological parameters and of the linear growth factor in the 
different models. On the other hand, as expected, the properties of collapsed 
structures do not appear to change significantly in Vector DE models, since 
no direct effect on the gravitational dynamics of particles is present in these 
cosmologies. Nevertheless, the growth rate of density perturbations shows a 
very peculiar shape in Vector DE cosmologies that clearly allows to distin- 
guish these models from standard Quintessence scenarios. 

The study of homogeneous DE models by means of their effects on nonlin- 
ear structure formation is presently entering the challenging era of precision 
cosmology, with a wealth of high-quality data expected for the near future. 
This implies the need to move from mainly qualitative assessments of the 
imprints of DE on the statistical and structural properties of self-gravitating 
systems to highly reliable quantitative estimations of the expected observa- 
tional footprints of each specific realization of a homogeneous DE field beyond 
A. Present and future simulations will then need to face the challenge of sig- 
nificanlty reducing statistical uncertainties mainly related to sample variance 
and to keep under control systematic effects due to numerical inaccuracies 
and most importantly to the yet poor understanding of sub-grid physical 
processes that are expected to heavily affect observable quantities at small 



scales (see e.g. Semboloni et al. , 2011). 



The former issue can be addressed by running larger cosmological simu- 
lations in terms of periodic box size, provided a sufficient mass resolution to 
resolve collapsed halos over a large enough range of masses can be achieved. 
Some attempts in this direction are presently being pursued with the Dark 



Energy Universe Simulations Series (DEUSS, Rasera et al. , 2010 Alimi et al 



2012) that aims to perform CDM-only simulations for the fiducial ACDM 
cosmology and for a few selected homogeneous DE models over simulated 
volumes comparable with the full observable Universe, employing a modified 



version of the AMR N-body code RAMSES (Teyssier, 2002). Such a chaL 
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lenge clearly requires highly sophisticated numerical tools with extremely 
high scalability and a dedicated pipeline for on-the-fiy data compression to 
maintain the volume of processed data still manageable. 

The latter issue, instead, does not show a similarly clear path towards 
possible solutions, and significant efforts will have to be made in the near 
future to refine our understanding of baryonic physics and astrophysical pro- 
cesses playing a substantial role in shaping the properties of cosmic structures 
at small scales, before these will be readily usable for cosmological studies. 



5. Simulating Dark Energy perturbations and interactions 

As it was discussed in Section [3} if the assumption of homogeneity of the 
DE field at sub-horizon scales is dropped, the effects of DE on the evolution 
of density perturbations and on the formation of linear and nonlinear struc- 
tures in the Universe are not confined anymore only in the Hubble friction 
term of Eq. |3| but can arise also through a direct contribution of the DE 
perturbations to the peculiar gravitational potentials experienced by mat- 



ter particles, as in Eq. 10, or even through additional interactions directly 



mediated by the inhomogeneous DE degree of freedom, as in Eqs. TTfT? 



Such scenarios require much more sophisticated algorithms to be properly 
implemented in N-body codes as compared to the simpler case of a homo- 
geneous DE field, which only requires to account for a modified expansion 
history through the correct Hubble function H{a). In the most general case, 
one should in fact devise algorithms capable to accurately solve a nonlinear 
Poisson equation like Eq. [13] for an arbitrary matter distribution with pe- 
riodic boundary conditions. This is an extremely challenging task, and the 
attempts to include such a sophisticated solver into N-body codes will be 
reviewed towards the end of this Section. However, this effort is in many 
cases not strictly necessary, as several specific DE models, although featur- 
ing sub- horizon perturbations and/or additional interactions, provide ways 
to directly relate the scalar field perturbations 6(f) to the matter distribution 
in the simulation box through simplified linear differential equations or even 
through algebraic relations. In this Section, I will review the results obtained 
with N-body simulations for these different classes of DE scenarios. 
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5.1. Non-universal couplings 

As no simulations have been performed so far for the case of non-interacting 
inhomogeneous DE cosmologies, as for instance the Clustering DE scenario 
introduced in Section Isl (see e.g. Sefusatti and Vernizzi, 2011), I will di- 



rectly move to review the results obtained in the last years for interacting 
DE models with non-universal couplings. These are scenarios for which ex- 
plicit screening mechanisms at small scales are not strictly necessary, which 
allows to significantly simplify the relation between the DE-mediated fifth- 
force and the matter distribution. A vast literature is available for a thorough 
description of the main features of this kind of interacting DE models, see 



cigalupi 


( 


2008 


); 



e.g. Amendola (2000 2004); Farrar and Peebles (2004); Pettorino and Bac 



(2009); Koyama et al. (2009); Baldi 



(2011b). As discussed above, for such models the scalar field perturbations 
can be simply related to the standard gravitational potential $ through an 
algebraic proportionality depending only on the coupling function. 



The first N-body simulations of Coupled DE models have been performed 



by Maccio et al. (2004) using a modified version of the AMR code ART for 
a range of cosmological models based on a DE scalar field with an inverse 
power-law potential of the form of Eq. |5] interacting with CDM only (i.e. with 
a vanishing coupling to baryons f3j, = 0) through a constant coupling function 
Pc in the range — 0.25. All the models were normalized to have the same 
amplitude of linear density perturbations at the present epoch, and evolved 
with a self-consistent background expansion history H{a). The early results 



of [Maccio et al.| showed that the fifth-force acting between CDM particles 
induces a bias between the amplitude of baryons and CDM perturbations, 
which is retained and amplified by nonlinear collapsed objects that show a 
reduced baryon content as compared to the standard ACDM case (see Fig. |4| 
left panel), and that the HMF at z = in Coupled DE models is practically 
indistinguishable from the ACDM case for a common ag normalization at 



the present epoch and can be accurately fit by the standard Jenkins et al. 



(2001) fitting formula. 



Another significant result of the early work of Maccio et al. is the dra- 
matic impact that the fifth-force was found to have on the inner slope of 
the halo density profiles - and consequently on the normalization of the 
concentration-mass relation - for the halos identified in the sample of their 
simulations: Coupled DE models were in fact found to produce highly over- 
concentrated halos as compared to ACDM, with a density profile approaching 
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Figure 4: Plots from Maccid et al. {2004) - Left: The ratio B{< R) between the baryonic 



and CDM overdensity enclosed in a sphere of radius R around the center of a massive halo 
extracted from the N-body simulations of different Coupled DE cosmologies, as a function 
of radius R. The linear bias between baryons and CDM that is visible at large distances 
from the center, due to the different gravitational dynamics of these two components, is 
significantly enhanced by nonlinearities in the inner part of the halo. Right: The density 
profile of a selected CDM halo forming in the standard ACDM cosmology (dot-dashed) 
and in various Coupled DE models with different coupling strengths. The effect of the DE- 
CDM interaction on the density profile appears dramatic in this work, with the maximum 
coupling /? = 0.25 giving rise to a steep power-law behavior of the profile with a logarithmic 
slope of ~ —2.3. Such trend has not been confirmed by the subsequent works of |Baldi 
et al.| ( |2010[ ) and |Li and Barro^ ( |2011[ ) (see Fig. [s] below) . 
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a power-law (and therefore not fit anymore by an NFW shape) with an inner 
logarithmic slope as low as ~ —2.3 for the largest coupling value (3c = 0.25 
considered in their work (see Fig. |4| right panel). Such result, which would 
have determined extremely tight constraints on the DE-CDM coupling as it 
would significantly worsen the cusp-core tension existing between numerical 
predictions and observations for the standard ACDM cosmology, was however 
not confirmed by later independent studies. 

In particular, the first adiabatic hydrodynamical simulations of Coupled 
DE cosmologies by Baldi et al. (2010) - performed with a modified version 
of GADGET - found essentially the opposite result for the same set of cos- 
mological scenarios: a mild reduction of the inner overdensity of halos for 
increasing values of the DE-CDM coupling jSc (see Fig. |5| left panel), with a 
consequent systematic shift of the normalization of the concentration-mass 
relation towards lower concentrations in Coupled DE as compared to ACDM. 
Baldi et al.| investigated further this issue by studying the impact of each 



individual modification of the standard ACDM dynamics implemented in 
their code, by running test simulations where each of these specific terms 
was artificially suppressed (see also Baldi, 2011a, for a systematic study of 
the different dynamical effects in interacting DE scenarios). As a result of 
this analysis, they concluded that the reduction of halo concentration was 
primarily determined by the effect of the friction term defined in Eq. 17 on 



the local particles dynamics: for a positive coupling (3^ > and a positive 
scalar filed velocity (p (which is what is realized for an inverse power-law run- 
away potential as the one assumed both by Maccio et al. and Baldi et al. ) 
the friction term Pc4'V acts as an effective drag (i.e. an "anti-friction") accel- 
erating particles along the direction of their motion. This corresponds to an 
injection of kinetic energy in virialized collapsed systems promoting the mi- 
gration of particles from inner to outer orbits, thereby adiabatically changing 
the virial equilibrium of the system towards more extended configurations of 
the halo core. This general result was then confirmed some time later by 
the independent collisionless simulations of Li and Barrow (2011) performed 
with a modified version of the AMR code MLAPM, that found a comparable 
shallowing of the inner density profile of CDM halos in Coupled DE models 
as the earlier results of IBaldi et al.l 

Besides the impact on the inner structure and concentration of collapsed 
halos, the work of Baldi et al. also investigated the specific baryon fraction 
in massive structures, finding that halos in Coupled DE cosmologies tend to 
have a significantly lower baryonic content than their ACDM counterparts 



33 




Figure 5: Plots from Baldi et al. (2010) - Left: The density profile of baryons (dot-dashed) 



and CDM (solid) for a massive halo extracted from the N-body simulations of the same 
Coupled DE models previously investigated by Maccio et al. (2004) (see Fig. |4] above). 
The trend of the inner overdensity of the halos as a function of the coupling j3c found 
in this work is the opposite of what previously claimed, with Coupled DE models giving 
rise to a reduction of halo concentrations. Right: The individual (open diamonds) and 
average (solid lines) halo baryon fraction in units of the cosmological baryon fraction for a 
sample including the 200 most massive halos detected in cosmological N-body simulations 
of the standard ACDM model and of various Coupled DE scenarios. Coupled DE induces 
a significant reduction of the halo baryon fraction at all masses. 
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(see Fig, [sj right panel), in good agreement with previous results. Finally, 
Baldi et al. studied in detail the evolution with redshift of the HMF, showing 



that both the analytical expression of Sheth and Tormen (1999) and the 



standard fitting function of Jenkins et al. (2001) reproduce with reasonable 
accuracy the HMF of Coupled DE cosmologies up to z ~ 2.5, provided the 
correct growth factor D^i^z) of each specific model is used for computing the 
theoretical halo abundance. 

The effects of Coupled DE models with a constant coupling to CDM 
on the high- 2; intergalactic medium, and in particular on the transmitted 
Lyman-a flux, has then been studied soon after with a series of radiative 



hydrodynamical N-body simulations by Baldi and Viel (2010), allowing to 



place new independent constraints on the coupling value of about ^ 0.15 
at 2(7 confidence level, while the impact on the correlation between CDM and 



galaxy distributions in clusters has been discussed in Baldi et al. (2011a). 



A related class of fifth-force models, where however the fifth-force is not 
necessarily associated with a DE degree of freedom but is rather assumed 
as a general additional interaction between massive particles, has been in- 



vestigated by Nusser et al. (2005), who ran cosmological N-body simulations 



including an additional fifth-force between massive particles, with the further 
complication that such fifth-force is assumed to be screened at large distances 
by a Yukawa suppression factor of the form exp(— r/r,,) in the fifth- force 
potential, with being a characteristic length scale defining the range of 
propagation of the scalar fifth-force 
work. 



see 



Gubser and Peebles 2004). In their 



Nusser et al. 



focused mainly on the effects of the fifth-force on the 
evolution of cosmic voids, finding that these specific fifth-force scenarios pro- 
duce a lower CDM and baryon density in voids as compared to ACDM (see 
Fig. [6| left panel), which is an appealing feature to address the longstanding 
problem of dwarf and irregular galaxies within voids being observationally 



too rare (Peebles, 2001). Similar studies have been subsequently performed 



also by e.g. Hellwing et al. (2010); Keselman et al. (2010). 



The same class of scenarios has then been tested also through N-body sim- 
ulations of individual galactic-size halos, focusing on the dynamics of dwarf 
satellites and on the effects of the fifth-force on the details of their tidal rem- 



nants. The first work of this kind, performed by [Kesden and KamionkowsE 
(2006), found very significant effects of the screened scalar fifth-force on the 



relative abundance of stars living in the leading and trailing tidal streams 
of gravitationally stripped dwarf satellites directly comparable to the Sagit- 
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Figure 6: Left, plot from Nusser et al. 



12005 ): the CDM and baryon distribution in a 
cosmological box of 50 Mpc/h aside, with and without an additional attractive scalar 
force with a Yukawa long-range suppression. The right upper corner shows the case of 
the standard ACDM cosmology, while the other plots display the particle distribution in 
the presence of a screened fifth-force. The most evident effect appearing from the plots 
is that cosmic voids are emptier than in the standard cosmological model. Right, plot 
from Kesden and Kamionkowski (2006): The stellar (red dots) and CDM (cyan dots) 
components of the remnant streams of a dwarf satellite tidally disrupted by its motion in 
the gravitational potential of a Milky Way-like spiral galaxy, as extracted from N-body 
simulations with and without (left lower plot) a scalar fifth-force. In the presence of a 
fifth-force the leading and trailing streams do no longer appear symmetric in their stellar 
content, which allows to place constraints on the strength of the fifth-force. 
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tarius stream. Such effects allowed Kesden and Kamionkowski to put very 
tight constraints on the maximum allowed value of the scalar interaction. 
However, a subsequent work by Keselman et al. (2009) found significantly 



different results for different choices of the initial conditions of the system, 
allowing for significantly larger values of the coupling without conflicting 
with direct observations of dwarf satellites tidal streams, although a further 



ditions chosen by Keselman et al. 



follow-up paper by Kesden (2009) challenged in turn the specific initial con- 



The possibility of a time-dependent coupling between DE and CDM, 
representing a more general class of interacting DE cosmologies than the 
constant coupling models simulated in the early works just discussed, has 
been included in N-body simulations for the first time in the work by [Baldi 
(2011b), that performed a series of adiabatic hydrodynamical simulations for 



different coupling functions including phenomenological parameterizations 
as e.g. f3c{a) oc a^^ and dynamical evolutions of the coupling such as e.g. 
/9c(a) oc exp [/3i0(a)]. In this work, also assuming a common normalization 
of the different models to the same cxs at z = 0, all the main basic analysis 
already performed for constant coupling models have been repeated, inter- 
estingly showing that the time variation of the interaction induces a whole 
range of new effects on structure formation processes that are in general ab- 
sent for the simplified case of a constant coupling. In particular, both the 
small-scale nonlinear power and the average halo concentrations, which can 
only be reduced as compared to ACDM within constant-coupling models, 
can instead show both trends - i.e. be either reduced or increased - for vari- 
able coupling scenarios, depending on the specific evolution of the coupling 
function (see Fig. [?]). This is due to the fact that besides the friction term 
(which as discussed above alters the virial equilibrium of collapsed objects 
by forcing them to expand) also the time evolution of the effective gravita- 
tional constant can modify the virial state of halos, and in particular for a 
coupling that grows in time this has the effect of favoring more concentrated 
configurations, thereby counteracting and possibly overcoming the opposite 
effect of the friction term. 



The effect of the linear amplitude normalization in interacting DE mod- 
els has been studied in a series of works using different assumptions for the 
simulations initial conditions. In particular, Baldi and Pettorino (2011) and 



Baldi (2012a) investigated the effects of Coupled DE models, both with con- 
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stant and variable coupling functions, on the expected number of massive 
clusters as a function of redshift, assuming a common normalization of the 
linear perturbations amplitude in the different models at high- 2; rather than 
at z = 0. These two works have shown how Coupled DE models consistent 
with the CMB normalization of the amplitude of density perturbations at 
last scattering systematically predict a larger abundance of massive clusters 
at high redshifts, as a consequence of the additional fifth-force that enhances 
structure formation thereby inducing a higher as normalization aX z = for 
cosmologies that start from the standard normalization at the redshift of 
last scattering {zis ~ 1100). Such result seems appealing to explain possible 
detections of extremely massive clusters at high- 2; that might result difficult 
to accommodate in the context of the standard ACDM scenario 



[see e.g. 



eFaTj [2009| |Rosati et aLj [2009t [Mortonson et al.\ [20TT1 



2011 



for an overview on this topic). In particular, Baldi (2012a) performed 



Jee 



Waizmann et al. 



the first N-body simulations of a specific type of Coupled DE models called 
"Bouncing" coupled DE, characterized by a constant coupling /3c and by a 
SUGRA self-interaction potential for the scalar field 0, resulting in a partic- 
ular dynamical evolution of the field that allows to match at the same time 
the normalization of linear density perturbations both at the last scattering 
surface zig ~ 1100 and at the present time, still allowing for significant de- 
viations from the ACDM behavior at intermediate redshifts. This work has 
shown that Coupled DE models of the "Bouncing" type allow to produce a 
significant excess of massive clusters at high redshifts without overpredicting 
the cluster counts in the local Universe, contrary to what can be achieved 
with standard Coupled DE models with constant or variable coupling and 
even with completely different approaches - such as e.g. primordial non- 
gaussianity - that have been invoked as a possible explanation for unexpect- 



LoVerde and Smith, 2011). 



Grossi et al. 


2007 


LoVerde et al. 


2008 



The different types of Coupled DE scenarios (constant /3, variable /3, 
"Bouncing") that have been studied through different (and often not mu- 
tually comparable) N-body simulations in the last years, have now been 
included in a large, systematic, and self-consistent simulations project with 
the aim to provide accurate and statistically significant numerical data for 
Coupled DE cosmologies to be readily compared with each other and with 
standard ACDM predictions. Such initiative goes under the name of the 
"Codecs Project" ( |Baldi| |2012cD and includes N-body and adiabatic hy- 
drodynamical simulations of a variety of cosmological scenarios, all sharing 



38 



10000.0 
1000.0 
100.0 



; 1 






ACDM = 




EXPQ1Qa2 - 




EXP015a3 - 




EXP010e2 ^ 




EXPQ10e3 : 




EXP015e3 


\ z=0.00000 

1 






Density profiles for Group nr. 



0.1 



1.0 10.0 
k [ h Mpc ' 1 



100.0 




ACDM 
EXP010a2 

EXP015a3 
EXP010e2 
EXPOIOeS 

EXP015e3 



z- 0.00000 

Mjojt'iCDM) = 3.B37S4e+14 h ' 



r[h ' Kpc] 



Figure 7: Plots taken from Baldi \201 1 b ) - The nonlinear matter power spectrum (Left) 
and the halo density profile (Right) as computed from a series of simulations of Coupled 
DE models with time-dependent coupling functions. The time evolution of the coupling 
induces additional effects on the evolution of structures and can modify the standard 
ACDM results (black lines) in opposite ways, depending on the specific time evolution of 
the coupling strength. 



the same WMAP7 (Komatsu et al. , 2011) cosmological parameters at z = 
(except for ag) and with a common normahzation of the hnear density per- 
turbations amphtude at the redshift of the last scattering surface. The post- 
processed numerical data of the "CoDECS Project" (such as nonlinear matter 
power spectra, halo and sub-halo catalogs, etc..) are made publicly available 
through a dedicated web-databas^and have already been used for a number 
of studies aimed at testing Coupled DE models against present or future ob- 



servational data. In particular, Lee and Baldi (2011 ) used these data to inves- 



tigate the impact of the DE-CDM interaction on the pairwise infall velocity 
of colliding galaxy clusters morphologically and dynamically comparable to 
the "Bullet" cluster (iMarkevitch et al.l 120021), finding that Coupled DE cos- 



mologies very significantly enhance the probability of high- velocity collisions. 



MaruUi et al. (2011) made use of the CoDECS public data to investigate the 



peculiar features of interacting DE in the clustering and redshift-space dis- 



tortions patterns of galaxies, while Beynon et al. (2012) computed forecasts 



for the weak lensing constraining power of the Dark Energy Survey (DES) 
and the Euclid satellite mission on the DE-CDM interaction, and ICui et al.l 



^http: / / www.marcobaldi.it / web / CoDECS 
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(2012) employed the same data to test the universahty of the HMF, finding 
evidence for deviations from the universal behavior at the level of about 10%. 



A completely different type of models belonging to the class of inhomo- 
geneous and interacting DE cosmologies with a non-universal coupling is 



given by the "Growing Neutrino" scenario (proposed by Amendola et al 
as a possible solution to the DE coincidence problem) , characerized by 



2008 



a constant and very large coupling to massive neutrinos {\/3i^\ > 50) while 
the other matter fields remain uncoupled. The evolution of perturbations in 
the Growing Neutrino model is characterized by a very fast growth of neu- 
trino density fiuctuations soon after the transition from the relativistic to the 
non-relativistic regime, which for realistic choices of the model's parameters 
happens around z ~ 4 — 6 (Mota et al. , 2008). As neutrino perturbations 



very quickly grow nonlinear, a full N-body treatment is required in order 
to properly follow the evolution of neutrino structures and of the relative 
gravitational potential at large scales. 

The first N-body simulations of the Growing Neutrino scenario have been 
performed by Baldi et al. (2011b) with a suitably modified version of GAD- 
GET for a model with coupling = —52 and a neutrino mass at ^ = of 
rriijQ = 2.4 eV. These early simulations have allowed to follow the formation 
of nonlinear neutrino halos at the expected scales of ~ 10 Mpc/h and larger 
down to 2; ~ 1, and to compute the backreaction effect of the gravitational 
potential associated with such neutrino lumps on the CDM distribution, 
showing a clear enhancement of the CDM bulk flow and an excess of CDM 
power at the largest scales available in the simulation box (^ 320 Mpc/h 
aside). The limitations of the newtonian approximation, which is generically 
assumed in N-body solvers, did not allow to run these simulations down 
to ^ = as the strong acceleration experienced by neutrino particles made 
the neutrinos relativistic again, with velocities comparable and eventually 
exceeding the speed of light at redshifts lower than 1. A signiflcant im- 
provement in this respect has been made with the completely independent 
numerical implementation developed by Ayaita et al. (2012) which explicitly 
includes a fully relativistic treatment of neutrino velocities ensuring that the 
speed of light limit for particles velocities be automatically fulflUed in the 
simulation. Additionally, the implementation of Ayaita et al. includes the 
effects of the local variation of neutrino masses and the backreaction of the 
formation of neutrino structures on the cosmic background expansion rate 
that were discarded in the earlier work of Baldi et al. (2011b). 
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5.2. Universal couplings and screening mechanisms 

I now move to consider the case of interacting DE models with univer- 
sal couplings, i.e. cosmological scenarios where the interaction between an 
inhomogeneous scalar degree of freedom that can be associated with DE in- 
volves all massive particles in the Universe. As discussed above, for such 
models an explicit screening mechanism capable to suppress the fifth-force in 
the local neighborhood of the solar system is generally required in order to 
avoid conflicts with local tests of General Relativity. Nevertheless, as a first 
order approximation for models with a coupling strength much smaller than 
gravity (i.e ^ 1) and with a sufficiently flat self-interaction potential, the 
issue of local recovery of standard gravity can be disregarded when focusing 
on structure formation processes at scales signiflcantly larger than the solar 
system itself. This is the case of scalar-tensor theories as e.g. Extended 
Quintessence models 



(see e.g. Perrotta et al. 


2000 




Baccigalupi et al. 


Pettorino and Baccigalupi 


201 IS 


) where a cosmic 



2000 



fleld playing the role of DE directly couples to gravity in the Jordan frame, 
corresponding to a universal coupling to matter in the Einstein frame (see 
again Pettorino and Baccigalupi 2008). Similarly to the case of Coupled DE 



models, also for Extended Quintessence theories it is then possible to sim- 
plify Eq. 13 and directly relate the strength of the extra flfth-force acting (in 



this case) among all massive particles to the standard gravitational potential 
through an algebraic equation relating the effective gravitational attraction 
experienced by massive particles to the standard Newton's constant. How- 
ever, differently from the case of Coupled DE, in Extended Quintessence 
models such relation directly depends on the sign of the effective coupling, 
such that the total interaction between particles can result both stronger or 
weaker than standard gravity, respectively for positive and negative values 
of the coupling. 



The flrst N-body simulations of Extended Quintessence scenarios have 
been presented in De Boni et al. (2011). They made use of the above- 



mentioned approximation of the effective flfth-force in terms of the standard 
gravitational potential for their modifled version of GADGET used to perform 
a series of radiative hydrodynamical cosmological simulations with gas cool- 
ing and star formation for a range of DE scenarios including also Extended 
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Quintessence models with both positive and negative couphngs. In their sim- 
ulations, |De Boni et al. always assumed a common normalization of linear 



density perturbations at last-scattering for all the different cosmologies, and 
focused on the hydrodynamical properties of massive halos corresponding 
to galaxy clusters. As a main result, |De Boni et al. showed that baryonic 



physics does not appear to be significantly affected by the additional in- 
teraction although the formation history of clusters and consequently their 
structural properties as well as their past record of star formation are altered 
by the DE phenomenology. Interestingly, they also found that both the stel- 
lar and gas content of relaxed, massive clusters is not significantly modified 
in cosmologies where a universal scalar interaction besides gravity is present, 
as compared to their ACDM counterparts. Such result provides a clear ob- 
servational way to discriminate between Extended Quintessence scenarios 
and models with non-universal couplings as the Coupled DE cosmologies 
discussed above, since for the latter the overall gas content of massive ha- 
los significantly chang GS clS db function of time as compared to the standard 
ACDM case (see Fig. g)- 

The accuracy of the approximation relating the extra fifth-force of Ex- 
tended Quintessence models to the standard gravitational potential has been 
explicitly tested by Li et al. (2011a) making use of a more sophisticated algo- 
rithm developed for modified gravity models that feature an explicit screen- 
ing mechanism (see below) implemented in a modified version of the AMR 
code MLAPM. Such algorithm is capable to solve the full nonlinear Pois- 
son equation [13] without resorting on any approximation for a wide range of 
functions F{6(j)) - including the case of scalar-tensor theories like Extended 
Quintessence - through a mesh-based iterative relaxation scheme. There- 



fore, in their simulations Li et al. could explicitly solve for the full scalar 
field perturbations 6(j){t, x) in a cosmological simulation box with periodic 
boundary conditions and derive the exact fifth-force acting on each particle 
within the simulated volume. By directly comparing the exact fifth-force 
computed with this algorithm to the one obtained by scaling the standard 



gravitational potential with the approximated relation adopted in De Boni 



et al. Li et al. showed that such approximation is highly accurate even 



for significantly larger coupling values than the ones investigated by De Boni 



et al. An explicit solution of the full nonlinear Poisson equation 13 is therefore 



not necessary for Extended Quintessence models. With their simulations Li 



let al.l also showed that several different observables like the nonlinear matter 
power spectrum, the halo mass function, and the concentration-mass rela- 
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tion are modified in opposite ways as compared to the standard ACDM case 
depending on the sign of the couphng. This is due to the fact that Extended 
Quintessence models are characterized by the superposition of two different 
effects related, respectively, to the modified expansion history and to the ex- 
tra fifth-force that characterize these cosmologies. In particular, the former 
effect tends to slow down the growth of linear density perturbations due to 
a faster expansion rate, while the latter can either enhance or suppress the 
growth of linear and nonlinear structures due to the larger or smaller effective 
gravitational constant. As a result, these two effects can either partially bal- 
ance each other (for the case of a positive coupling, i.e. an enhanced effective 
gravity) or conspire towards a significantly slower growth of perturbations 
(for the case of a negative coupling). In the former case, linear perturbations 
are only mildly suppressed or almost unaffected, while in the latter they re- 
sult significantly suppressed. At nonlinear scales, however, the time variation 
of the extra fifth-force becomes the dominant effect giving rise to an excess 
of power and a significant increase of the concentration of CDM halos for 
models with a negative coupling, as these feature a positive derivative of the 
effective gravitational constant at low while for positive couplings (i.e. a 
decreasing effective gravitational constant) the opposite effect arises. For all 
cases, the HMF is found to be suppressed as compared to ACDM. 

Once taking into account the fact that |Li et al. adopted the same initial 
conditions for all their different DE simulations, thereby implicitly imposing a 
common normalization of the linear perturbations amplitude of all the models 
at some intermediate redshift Zi between last scattering and the present time, 
their results appear to be in general qualitative agreement with the earlier 
outcomes of De Boni et al. (2011). 



As a follow-up to their early work, De Boni et al. (2012 ) extended the anal- 
ysis of their simulations to a detailed investigation of the concentration-mass 
relation for clusters in DE models including Extended Quintessence scenarios, 
finding again that the sign of the time derivative of the effective gravitational 
constant drives the shift in the normalization of the concentration-mass re- 
lation, consistently with the outcomes of Li et al. (2011a). This effect is 
somewhat similar to the one detected for Coupled DE models with a vari- 
able coupling (Baldi, 2011b, see above and Fig.js]) where the time variation of 
the effective gravitational constant alters the virial equilibrium of collapsed 
objects inducing a contraction or an expansion of the halos and consequently 
an increase or a decrease of their concentration parameter as compared to 
ACDM. 
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The most general case of a universal interaction between a cosmic inhomo- 
geneous scalar and matter fields in the Universe corresponds to cosmological 
models where nonlinearities in the function F{6(j)) appearing in Eq. 13 induce 
large spatial fluctuations in the scalar field, capable to provide an efficient 
screening of the extra fifth-force at small scales even for effective coupling 
values of order unity, i.e. for a strength of the fifth-force comparable to grav- 
ity. If significant nonlinearities in the function F are present, in fact, it is no 
longer possible to discard the term F{6(j)) in Eq. 13 and approximately relate 
the scalar field perturbations 50 to the matter perturbaitons SpM through a 
standard linear Poisson equation. This is the case of modified gravity mod- 
els as e.g. f{R) theories, Symmetron fields, or higher-dimensional theories 
of gravity that were introduced in Section [3] above. In these classes of mod- 
els, then, it is strictly necessary to solve the full Eq. [13] in order to compute 
the actual fifth-force acting on massive particles at different positions, with- 
out resorting on any further approximation, as the behavior of the extra 
force will be different in different environments due to the explicit screening 
mechanisms defined by the details of the function F{5(j)) and of the coupling 
function This is an extremely challenging task in itself, which requires 

dedicated algorithms to be included and interfaced with standard N-body 
solvers, and that sensibly increases the computational cost of large N-body 
runs for these scenarios. 

The field of cosmological simulations of modified gravity models, and in 
general of scalar field cosmologies with explicit screening mechanisms, is still 
rather young, although in recent years the efforts to develop competitive and 
versatile N-body codes for this class of scenarios have been remarkable. A 
proper review of such field is then probably still premature, since it is only 
very recently that independent simulation codes have started to produce 
broadly consistent results for some specific realizations of modified gravity 
theories, and a proper comparison of different algorithms for cross-checking 
and mutual validation has not yet been performed. Nevertheless, the amount 
of work invested by several different research groups in developing and testing 
such implementations over the last few years is definitely worth a mention, 
along with the clear prediction that in a relatively short timescale a large 
amount of robust and significant results in this field will be achieved through 
a systematic program of numerical investigations. I will therefore provide 
here only a very brief and general overview of this field, and leave to a future 
time a more thorough discussion. 
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Figure 8: The effect of a series of f{R) modified gravity tlieories on the nonhnear matter 
power spectrum as compared to the standard ACDM cosmology within General Relativity. 



The plots are extracted from the first simulations of f{R) models by jOyaizu et al. (2008) 
(Left) and from the more recent results of Li et al. (2012b I (Right). The enhancement in 



the nonlinear matter power spectrum due to the additional fifth-force related to the mod- 
ified gravity theory is significantly suppressed at small scales by the Chameleon screening 
mechanism, in a way that appears consistent between these two different implementations 
of f{R) gravity models in N-body simulations. 



The first attempt to run N-body simulations for modified gravity in tlie 



form of f{R) tlieories was made by Oyaizu (2008); Oyaizu et al. (2008); 



Schmidt et al. (2009), who made use of an iterative relaxation scheme on a 



fixed cartesian grid within a mesh-based N-body code to solve Eq. [13] and 
self-consistently compute the motion of particles in the presence of a modified 
gravity law given by the superposition of the standard gravitational interac- 
tion and a screened fifth-force. Their results showed for the first time how 
the screening mechanism strongly suppresses the effects of the fifth-force at 
small scales, and in particular in the inner parts of collapsed halos where 
the suppression is most effective (see Fig, [sj left p anel) . A similar numerical 
approach was followed soon after by Li and Zhao (2009); Zhao et al. (2010); 



Zhao et al. (2011 ) who implemented a Newton-Gauss-Seidl iterative solver for 



the scalar field nonlinear Poisson equation on the adaptive grid of the AMR 
code MLAPM, allowing for a significant improvement of the code resolution 
as compared to the earlier implementation of Oyaizu in regions where the 
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screening mechanism is most efficient. The relatively small simulations that 
could be carried out with such AMR implementation of modified gravity gave 
results in good agreement with previous findings, and triggered a significant 
number of post-processing analysis aimed at investigating possible charac- 
teristic features of modified gravity models in various observables. These 
include e.g. the statistics of CDM halos and voids (as in Zhao et al. , 2010 



Zhao et al. 2011; Winther et al. 2011 Li et al. 2011c), the geometrical and 



dynamical properties of virialized halos (Lombriser et al. 2012 Lee et al. 



2012 Llinares and Mota 



Li et al. 



2012 ) or the apparent variation of the fine-structure 



2011b) as well as a number of applications to dif- 



constant (as m 

ferent types of screening mechanisms such as the Dilaton (Brax et al. , 2011) 
and the Symmetron (Davis et al. , 2012 ). The same implementation has been 
recently ported by Li et al. (2012b) into the hydrodynamical AMR N-body 
code RAMSES and has been named ECOSMOG standing for Efficient COde 
for Simulating MOdified Gravity. Such code overcomes several shortcomings 
of the previous MLAPM implementation concerning the parallelization strat- 
egy and the multi-grid refinement for solving the nonlinear Poisson equation, 
which makes the ECOSMOG code more suitable for large simulations with 
high mass resolution. The first series of simulations performed with such code 
have been focused on f{R) models characterized by a Chameleon screening 
mechanism as well as on different types of screening as e.g. Symmetron- and 
Dilaton-type modified gravity theories (see e.g. Brax et al. , 2012), and are 
now starting to be post-processed and analyzed with a particular focus on 
the effects of modified gravity theories on the nonlinear matter and velocity 



power spectra (Li et al. , 2012a) and their connection with direct observables 



such as the detailed pattern of redshift-space distortions in wide galaxy sur- 
veys (Jennings et al. , 2012), finding consistent results with previous works 



(see Fig. [8| right panel) . 

The number of independent implementations of modified gravity theories 
into cosmological N-body codes and the range of accurate predictions for di- 
rect observable quantities that are being produced with such last generation 
of N-body solvers seems encouraging in view of the needs of large upcom- 
ing surveys aimed at investigating the nature of the dark sector and of the 
gravitational interaction at cosmological scales. 
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6. Simulating large-void cosmological models 



As a last class of models of the cosmic acceleration beyond A that have 
been investigated by means of N-body simulations it is worth to mention the 
case of large- void inhomogeneous cosmologies introduced in Section [3j The 
only cosmological simulations of such models so far have been performed by 



Alonso et al. (2010) using the TreePM code GADGET. Such simulations did 



not require any specific modification of the public version of the GADGET-2 
code but rather a modification of the initial conditions generator to include 
the gravitational potential of a large void in the computation of the initial 
particles' displacements. This has been done by suitably modifying the code 
2LPT based on 2nd order Lagrangian perturbation theory to include in the 
initial conditions a series of large voids with different values of the structural 
parameters and Ar/r^, that have then been evolved using the standard 
newtonian N-body approach. One of the most relevant results of such inves- 
tigation has been to demonstrate that standard N-body codes are suitable to 
correctly follow the nonlinear evolution of the density void (i.e. to correctly 
predict its evolution up to a density contrast 5 ~ 1 at the present time) 
without requiring any general relativistic modifications. Also, the same work 
showed that the linear density contrast in such inhomogeneous cosmologies 
is always very close to the standard ACDM prediction. Such results sup- 
port the employment of large cosmological newtonian N-body simulations 
also to investigate scenarios for which the basic assumption of large-scale 
homogeneity encoded by the Copernican principle does no longer hold. 



7. Conclusions 

The last decades of investigation in the fields of Cosmology, Astrophysics, 
and Particle Physics, have provided us with a clear and undeniable quantifi- 
cation of our ignorance: about 96% of the energy density in the Universe is 
made of particles and fields that keep eluding all our efforts of detection and 
identification. In such context. Dark Energy and Dark Matter are therefore 
simply labels that allow us to organize and classify the limited observational 
knowledge that we have been growing through the years about such "dark 
side" of the Universe. Although the simplest and most widely accepted cos- 
mological model - that associates Dark Energy to a cosmological constant 
and Dark Matter to a family of Weakly Interacting Massive Particles beyond 
the standard model of particle physics - is presently consistent with all our 
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available observational data, its theoretical roots are difficult to accommo- 
date in the context of General Relativity and Quantum Field Theory, and 
alternative scenarios keep being proposed almost on a daily basis since more 
than a decade. 

We are therefore accumulating an increasing number of alternative cos- 
mological models that aim at providing a solution to the mystery of the fun- 
damental nature of the dark Universe, which are often barely distinguishable 
from each other in their predictions concerning the background evolution of 
the Universe or the growth of linear density perturbations. Trying to exploit 
also the nonlinear regime of structure formation as a possible way to discrim- 
inate among different cosmological scenarios and as a source of observational 
information about the nature of the dark Universe is therefore becoming a 
necessary further step in the connection between theory and observations 
in cosmology. Such a step however requires to make use of large numerical 
simulations as the nonlinearities involved in the problem prevent to drive 
any reliable conclusion based only on analytical tools. This need has driven 
the wide range of efforts that have been made in the last years to develop, 
test, and ultimately apply new and highly sophisticated algorithms within 
N-body codes to self- consistently simulate the evolution of cosmic structures 
in the context of different and competing cosmological scenarios. 

In this Review, I tried to provide a broad overview on the results of such 
new and rapidly developing research field, mainly focusing on cosmologi- 
cal N-body simulations of Dark Energy models alternative to the standard 
ACDM scenario. After briefly reviewing the history of the role played by N- 
body simulations in estabhshing the present standard cosmological model, I 
provided a broad (and necessarily incomplete) overview of the different Dark 
Energy scenarios that are presently being considered as possible competi- 
tors to the standard model. In doing so, I classified Dark Energy models 
in two different categories defined by the clustering properties of the Dark 
Energy field (whatever this field might be) at sub-horizon scales, dehber- 
ately avoiding any attempt to make a fundamental distinction between Dark 
Energy and Modified Gravity scenarios. In fact, no fundamental distinction 
between a Dark Energy component in the stress-energy tensor of the Uni- 
verse and a modification of the laws of gravity is possible when one allows 
for density perturbations and direct interactions of the Dark Energy field. 
The same classification, however, results also particularly useful to discuss 
the specific modifications that have to be implemented in cosmological N- 
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body algorithms to account for the characteristic features of different Dark 
Energy models. In fact, depending on the clustering properties of the Dark 
Energy field, structure formation can be affected either only through a mod- 
ified background expansion history, or also by additional forces related to 
the local Dark Energy density perturbations. The numerical investigation of 
these two distinct possibilities through N-body simulations has been the main 
focus of the present Review. A further possibility, which does not belong to 
any of these two main categories, is given by models that relate the observed 
accelerated expansion of the Universe to a local deviation form homogeneity. 
Such option requires a different kind of numerical implementation, and has 
been discussed separately. 

In the second part of this work, I attempted a general overview of the 
main investigations that have been performed in the last decade using N-body 
simulations of non-standard Dark Energy models, following the general clas- 
sification summarized above, and focusing the discussion only on the main 
outcomes of the various studies rather then on their technical details. Due 
to the complexity of the field, and to the wide range of different cosmolog- 
ical models, N-body codes, and normalization choices assumed in different 
studies, the description of most of the mentioned works has necessarily been 
incomplete and oversimplified, but I tried to highlight the most relevant re- 
sults obtained by different research groups and to which extent such results 
have been subsequently confirmed or disproved by other independent inves- 
tigations. In any case, I tried to provide an extensive list of references to 
address interested readers to the relevant literature. 

In this broad overview, I mainly focused on homogeneous Dark Energy 
models and on various types of interacting Dark Energy cosmologies, that 
are the classes of models for which a wide number of independent cosmo- 
logical simulations have been carried out so far, providing consistent results 
which can then be considered sufficiently robustly estabhshed. In the last 
part of the Review, however, I provided also a brief and necessarily incom- 
plete summary of the main efforts that have been put in place in the last 
years to develop N-body codes for various types of Modified Gravity models, 
that according to the above-mentioned classification scheme correspond to 
inhomogeneous Dark Energy models with a universal screened interaction to 
massive particles. The field of N-body simulations of Modified Gravity mod- 
els is presently very active and has shown an impressive progress in the last 
couple of years, but the very recent development and application of most of 
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the presently available codes together with the lack of a direct comparison of 
different implementations make a proper review of this specific field probably 
still premature. 

The number of different and independent efforts aimed at developing suit- 
able numerical tools to push the comparison between theoretical models of 
the dark side of the Universe and direct observations deep into the nonlin- 
ear regime of structure formation has been continuously growing in the last 
decade. Despite the difficulties related to the intrinsic nonlinear nature of the 
processes under investigation, the field of N-body simulations of Dark Energy 
and Modified Gravity models seems to be rapidly developing and promises to 
provide highly constraining predictions for a wide range of presently viable 
and competing cosmological scenarios. The additional complications, which 
have not been discussed in the present work, related to possible degenera- 
cies with other physical processes in place at the same nonlinear scales at 
which present N-body simulations of non-standard cosmologies are making 
their most valuable predictions, will necessarily have to be taken in full con- 
sideration in the future. Also, a synergy between a proper implementation 
of Dark Energy models and a better understanding of baryonic physics will 
be required to obtain the level of accuracy which is demanded by the next 
generation of observational surveys. 

Although an impressive range of results have been obtained in the last 
decade from N-body simulations of non-standard cosmological scenarios, we 
are only at the beginning of a long and challenging path for numerical cos- 
mology, and the present work is probably only the first of a long series of 
Reviews in this new and exciting field of research. 
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